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relevant weak system-bath couplings that is applicable even in the thermodynamic limit. We identify conditions 
under which thermalization happens and discuss the underlying physics. Based on these results, we also present 
a fully general quantum algorithm for preparing Gibbs states on a quantum computer with a certified runtime 
and error bound. This complements quantum Metropolis algorithms, which are expected to be efficient but have 
no known runtime estimates and only work for local Hamiltonians. 
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How do thermal quantum states - cornerstones of a descrip- 
tion in canonical ensembles in quantum statistical physics - 
arise from the underlying theory of quantum physics? This 
question, a long tradition as it obviously has, is in many ways 
still surprisingly wide open. Indeed, much progress was made 
only recently ifn- fTTl : this is motivated and triggered both 
by new mathematical Il5lj9l [TTl [T2l and numerical [ 1 3 1 tech- 
niques becoming available, as well as by new experiments 
with quantum many-body systems in nonequilibrium fl4l . 

In this work we present a set of precise sufficient condi- 
tions for the emergence of Gibbs states from the underlying 
microscopic theory of quantum mechanics. Our results go be- 
yond previous approaches in that they apply in a physically 
relevant weak coupling limit and constitute the key insight 
leading to the invention of a quantum algorithm that prepares 
Gibbs states with certified precision and runtime. 

The three ingredients that enter the standard textbook proof 
of the canonical ensemble in classical statistical physics are: 
(i) the equal a priori probability postulate (also known as mi- 
crocanonical ensemble) and an equilibration postulate (such 
as the second law), (ii) the assumption of weak coupling, and 
(iii) an assumption about the density of states of the bath, 
namely, that it grows faster than exponentially with the energy 
and that it can be locally well approximated by an exponen- 
tial [15J. Here each of these steps is translated to the pure 
state quantum statistical mechanics approach ITIjS). In par- 
ticular (i) can be replaced by either a typicality argument, or 
a statement about dynamical relaxation that follows directly 
from quantum mechanics and (ii) is made precise by proving 
a novel perturbation theorem that has applications far beyond 
the scope of the present article. 

Our new technical results allow us to design a quantum al- 
gorithm preparing Gibbs states with explicit error and run- 
time bounds, invoking a new variant of phase estimation. Our 
algorithm complements another algorithm with certified run- 
time that was proposed in Ref. ||T6| and recent developments 
on quantum Metropolis algorithms fTTI . 

Setting and notation. We consider a system S weakly 
coupled to an environment B. The Hilbert space reads T-L — 
Hs <E)'Hb, where Hs and "H^ are the Hilbert spaces of the 



subsystem and the "bath" (with finite dimensions ds and ds)- 
The evolution of the total system is governed by the Hamil- 
tonian H = Ha + V, with eigenvalues and eigenvectors 
{Ek} and {jE'fc)} consisting of an uncoupled Hamiltonian 
Hq = Hs + Hg, with eigenvalues and eigenvectors {E^j^^} 
and {{E^'^)}, and a coupling Hamiltonian V. We give con- 
ditions under which the reduced state ipf — Tr^ ^pt, with 
"0* — IV"*) {i^t\, of the subsystem S relaxes for most times to 
a Gibbs state Pcibbs •= c^^^^ / Tre^^^^ with inverse tem- 
perature P under unitary time evolution \ipt) = Q-'^Ht |^^^^ 
By this we mean that for most times their trace distance 
I?(i/'f , Pgjjjj^g), which measures the physical distinguishabil- 
ity 1211 . is small. Note that the decomposition of a given H 
into Hs, Hb, and V is not unique. This freedom can be 
used to optimize the bounds in our results, and the correct 
Hs naturally results from this optimization. We assume that 
the Hamiltonians H and Hq are non degenerate such that time 
averaging and dephasing in the eigenbasis give the same re- 
sult 
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Whenever an expectation value equilibrates, it does so to the 
expectation value in w flT] . 

"Natural thermalization" : Conditions for Gibbs states to 
appear In this section we go through points (i)-(iii). The 
final conclusion is summarized in Corollary [T] The central 
point of the argument is a novel perturbation theorem that re- 
lates spectral projectors of weakly interacting and noninteract- 
ing Hamiltonians in a physically relevant weak coupling limit. 
It allows us to connect results on dynamical equilibration and 
measure concentration with classical counting arguments and 
thereby prove a set of natural sufficient conditions for ther- 
malization in quantum mechanics. 

A stepping stone in the argument will be states that have 
an energy distribution that is flat in an interval [E, i? + A] 
and vanishes otherwise. We indicate such states, and their 
dephased states, by a subscript □ like in V'n or wn and call 
them rectangular states. This class of states includes both 
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mixed states (in particular the microcanonical state wn) and 
pure states and thus, because of the freedom to choose the 
phases, usually also initial states that can be locally out of 
equilibrium. 

The equal a priori probability postulate (i) can be replaced 
by a typicality argument using results from Refs. lH] El 13. 
In Ref. 1 1 1 powerful concentration of measure techniques are 
used to show that almost all states from a microcanonical 
subspace corresponding to a microcanonical energy window 
[E, E + A] locally look like the reduction of the correspond- 
ing microcanonical state; i.e., Wp) is small for all but 
exponentially few of the states ip from the subspace, where 
LOn is the microcanonical state on the subspace. Alternatively 
one can use the results concerning the dynamics of states with 
a high effective dimension of Refs. ||5]|6|. Under one assump- 
tion on the spectrum of the Hamiltonian (nondegenerate en- 
ergy gaps), it is shown that all reduced states on small sub- 
systems of such states tend to the time-averaged equilibrium 
state uj^ and stay close to it for most times. In many-body 
systems, natural initial states have a high effective dimension, 
and this is provably true for all but exponentially few states 
from a microcanonical subspace |5|. 

The delicate issue, which has up to now not been addressed 
in the literature in a general and rigorous way, is the weak cou- 
pling approximation (ii) |6, 22). The problem is that due to the 
exponential growth of the Hilbert space dimension and the at 
most polynomial growth of the energy content, the spectrum 
of the noninteracting Hamiltonian Ho becomes exponentially 
dense with increasing bath size. Therefore, the perturbative 
limit, in which the coupling V is weak compared to the gaps 
of the noninteracting Hamiltonian H^, and in which it can 
be guaranteed that the energy eigenvectors \Ek) of the full 
Hamiltonian H = Hq + V are close to product states, is ar- 
guably not the physically relevant weak coupling limit. Even 
worse, in this limit memory effects provably prevent thermal- 
ization ||8|. As in the classical setting, a coupling should be 
considered to be weak as long as it does not change the total 
energy in a noticeable way. That is to say, the energy stored 
in the interaction is much less than our (microcanonical) un- 
certainty about the energy of the system, i.e., ||y||oo ^ A, or 
for thermalizing systems much less than the thermal energy 
1//3. This is the relevant weak coupling limit in which we 
prove equilibration towards a Gibbs state. We do this by re- 
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of the dimensions of the supports of uip^'^ and Wn, and Q.^ is 
the total number of eigenstates of H and Hq in the intervals 

[E, E + e] and [E + A - £,E + A]. 

The theorem shows that, for any two initial (possibly 
pure) states that have a flat energy distribution in the interval 
[£', i? + A] with respect to the Hamiltonians Hq and H with 
||y||oo ^ A, the distance of their reduced dephased states 
Wp''^' and is small. In particular, assuming an approxi- 
mately constant density of states such that f7£/(2 fimax) ~ 
2e/A andMV(2fii„ax) ^ ||F||oo/A, the best choice for e is 
£ w v^||F||ooA/2 which gives 
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In cases with an exponential density of states, for which we 
will get equilibration towards PQib^g oc e~^^^, we can guar- 
antee that I?(ajp, cjp is small whenever ||V^||co ^ 
(compare Appendix H). 

Proof. First note that by monotonicity of the trace distance 
and the triangle inequality 
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where G and F are the projectors onto the sup- 
port of wn and Wp''' respectively and ilmin/max — 
min / max(rank(G), rank(F)), and Ail — fimax — ^^min- It 
remains to bound ||G-F||i.LetG = l- GandF = l- i^; 
then G - F = GF - GF, and thus ||G-F||i < ||GF||i + 
||GF||i. To bound ||GF|ji we decompose G — Gi + Ge into 
an interior part Gi which is the projector onto the eigenstates 
from the interval [E + e,E + A — e] and the exterior part 
Ge and find ||GF||i < |lG,F||i + ||Ge||i (seeFig.[l}. Using 
the inequality || • ||i < rank(-)|| • ||oo, submultiplicativity of 
the rank, and that rank(Gi) < ilmax this can be recast into 
||GF||i < r2max||G,F||oo + rank(Ge). Finally, from Theo- 
rem V.II.3.1 in Ref. HI it follows that ||G,F||oo < ||V"||oo/e- 
Repeating the argument for ||GF||i, introducing the notation 
fig = rank(Ge) + rank(Fe), and putting everything together 
gives the desired result. □ 

The level counting argument (iii) - with is ultimately the 



lating the dephased/microcanonical state ojn to the state cjn ^, 
dephased with respect to the non interacting Hamiltonian, for reason for the exponential form of Pcibbs 
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which we can easily perform the partial trace to obtain ajn'""* 
and thereby an approximation to Wp. 

Theorem 1 (interacting vs. noninteracting case). Letuj^'' and 
cjn be the dephased/microcanonical states belonging to the 
interval [E^ _E + A] with respect to Hq and H = Hq + V; 
then for every e < A/2 

P(.,^.n^("))<P(.n,a.("))<& + ^J±^, (1) 

£ Z i 'rnax 



ries over to the quantum case in a straightforward way in 
the absence of coupling between system and bath |2, 6|, 
and with a bit more work one can also obtain a rigorous 
trace norm error bound. If the number of states of the 
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where r!„,ax and AQ are the maximum, and the difference, ^'^'^ ^ ■ Definition of the projectors used in the proof of Theorem[T] 
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bath n^{E^) in the interval [E^,E^ + A] is such that 
the proportion fl^iE - E^)/J2i ^li^ - Ef) is close to 
Q-P^k I e~P^i for the given E and A and some (5, then 

the distance of 2^(^^n^°'', Pcibbs) small. This can be guar- 
anteed under a set of natural assumptions that are satisfied 
by a wide range of natural quantum many-body systems and 
that resemble the ones commonly used in classical statisti- 
cal physics, such as an exponential increase of the density of 
states (Appendix A). In particular, for a bath consisting of m 
noninteracting spin- 1/2 particles with a slightly varying on 
site field strength and average local energies of and rj one 
finds (Appendix B) 



l|gsll 



-l) + C 
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with C exponentially small in the bath size. We will later use 
this bath in our algorithm. In summary, Eq. Theorem [T] 
and the results on dynamical equilibration and random states 
from the unitary invariant measure derived in Refs. |lj|5J lead 
to the following conclusions: 

Corollary 1. (Kinematic) Almost all pure states ij) from a mi- 
crocanonical subspace corresponding to an energy interval 
[E^ + A] of a weakly interacting, sufficiently large quantum 
system are locally close to a Gibbs state in the sense that for 
every gaps(i7o) <C e < A/2 the probability that 
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drops of exponentially with f2„ 



(Dynamic) Moreover, if 



the Hamiltonian in addition has nondegenerate energy gaps 
all initial states '!/'n,t=o. even those locally out of equi- 
librium, with aflat energy distribution in the interval, locally 
equilibrate towards Pqi^^i^s in the sense that 



V(^s si'^)< 
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Both inequalities are robust against deviations from the rect- 
angular distribution. If the bath has an exponentially increas- 
ing density of states only a region of bounded variation fol- 
lowed by a sharp cutoff towards higher energies should be 
sufficient (for details see Appendix C). 

"Artificial thermalization" : A quantum algorithm for 
Gibbs state preparation. It follows from Eq. Q and The- 
orem[T]that all one has to do to prepare a Gibbs state is to pre- 
pare a state close to or ajp''' on a suitable combination of 
system plus bath. This is the central idea behind the quantum 
circuit shown in Fig.|2] which prepares thermal states without 
using any knowledge about the eigenstates of the Hamilto- 
nian. 

Quantum algorithms that prepare thermal states have sev- 
eral advantages over classical simulation methods: Quantum 
Monte-Carlo methods offer a way to, for example, estimate 



correlation functions of thermal states on a classical com- 
puter However, such methods are restricted to certain types of 
Hamiltonians as they suffer from the sign problem. A proce- 
dure that certifiably prepares Gibbs states in a quantum com- 
puter does not only overcome the sign problem, but moreover 
makes it possible to use thermal state in experiments address- 
ing questions of nonequilibrium dynamics in quantum simu- 
lators, for example to study quenches. 

Our algorithm requires two registers (see Fig. [2]i. The first 
register R consists of r qubits initially in |0) and is used to per- 
form quantum phase estimation. The second register Q holds 
the quantum system plus bath and is put into a rectangular 
state by performing the following steps: 

1. Initialization. The register Q is initialized into the com- 
pletely mixed state pi ^ \ Y1=i \Ek){Ek\ ® ■ 

2. Partial quantum phase estimation. A new form of quan- 
tum phase estimation is performed, which comprises three 
steps: the application of r Hadamard gates on the qubits of R, 
the application of r controlled-C/ operations (with U raised to 
successive powers of two), and an inverse Fourier transform 
on R. After this operation, the state of the registers is 



P2 



2^-1 d 

EE 

s,s'=0 fc = l 



as{vk)al,{ipk)\Ek){Ek\®\s){s'\, (7) 



where Lpk :- 

1 l ~cxp(27r i(2''yfc- 
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Note that 



and asi-^k) ■= 
\as{ifk)\'^ is a prob- 



2^ l-cxp(27ri(ipfc-s/2'-)) ' 

ability distribution that becomes more and more peaked 
around s/2'' as r increases. 

3. Measurement. Measuring the first q qubits of R, some 
binary string s* of length q is obtained and the system is left 
in the state 
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J2oisiVk)a:,i^k)\Ek){Ek\^\s){s'\, (8) 



=s,A, k=l 



where A* := 2'"^' is the number of states of the ancilla reg- 
ister R compatible with the measurement. By choosing r one 
can determine the width A — ||iJ||oo2~'^A* of the rectangu- 
lar state that is prepared. The measured value of determines 
the energy E — ||iJ||co2~'^s* of the rectangular state, and 
thereby the inverse temperature /3 of the Gibbs state. To ther- 
malize the subsystem at some particular temperature, the pre- 
vious steps must be repeated until the desired energy is mea- 
sured. The number of runs increases exponentially with the in- 
verse temperature (3. This prevents us from preparing thermal 
states at very low temperatures (see Appendix D). This is not 
a deficit of the algorithm, for otherwise QMA-hard problems 
(the quantum analog of NP lfT9]| ) could be efficiently solved. 
Any general algorithm will presumably have this feature |fT9l . 
The final state of Q is 



d /(s. + l)A. 

UJQC := Tr^pa cx ^ ^ \asiipk)\'^ ) \Ek){Ek\ 

fe=l \ s=s,A, 



(9) 
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Figure 2. Quantum circuit that generates a dephased rectangular state 
ojn''^ . I is tlie initialization gate, H are Hadamard gates, Ur = , 
U = exp(- i/Zo/ll-ffolloo) with Hq = Hs + Hb, and F'' is the 
inverse Fourier transform. 



For large enough r, this state is close to the desired state wn 
with^; = ||i?||oo2"«s* and A = ||i?||oo2"''A*. The precise 
deviation of wg^ from pgi^^^. 



I^KcPGibbJ <2?(^QC,w^°')+^?(w^^"^PGibbs)• (10) 

depends on the density of states of system plus bath. A good 
candidate for the bath is the system of m non interacting spin- 
1/2 particles discussed before (Appendix B) and we give ex- 
plicit results for the errors and the complexity for this bath: 

Algorithm. For any chosen A > 0, any given inverse temper- 
ature l3 and system Hamiltonian Hs, the algorithm presented 
in Fig. [2] using the bath with m spin-1/2 particles and en- 
ergy scale 1] = \/X/m\\Hs\\oo discussed before (Appendix 
B), prepares the system S of n qubits in a state within trace 
norm distance bounded by 



,■3(0) „S 



H^'qc. PGibbs) < 2^-''+^ (1 + ln(2'-^)/^2) 

,f+/3|lffsllo= + 



(11) 
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with C exponentially small in m, to a Gibbs state Pq^y^i^^ with 
a temperature in the interval [/? — 5/3, /3 + (5/3], where 



(5/3 < 2 



2-q 




(12) 



spins. In practice, in absence of an oracle for the Hamiltonian 
of the system, the error produced to perform the U gate carries 
a second source of error that comes from the Trotter-Suzuki 
approximation. Nevertheless, this eiTor can be suppressed at 
a polynomial cost for local Hamiltonians lfT6ll20l . 

The two contributions to the trace distance error ( [TO| i are 
computed in Appendices B and F, the average number of runs 
is computed in Appendix F and the error in the temperature 
comes from the discrete nature of the energy measurement via 
quantum phase estimation and is calculated in Appendix E. 
As is clear form Fig. [i] we need X]t=o ^ "^^ ^ 
gates, and the part of the circuit that does not correspond to 
the controlled time evolution, i.e., the initialization and the 
inverse Fourier transform, only requires the implementation 
of ?T. + TO + 2r Hadamard gates, r controlled single qubit gates, 
and n + m + q single qubit measurements. 

For a fixed large A, the error of the algorithm can be made 
small by choosing an ancilla register of size 0(/3^||ifs||^). 
The exponential scaling of the runtime with /3||i/5||oo caused 
by this is not a deficit of the algorithm, as any general efficient 
algorithm would contradict hardness results such as the local 
Hamiltonian problem |19|. Unlike our approach. Metropolis 
algorithms | HJ are expected to be efficient in some cases, but 
no rigorous runtime estimates are known and they are only 
applicable to local Hamiltonians. An interesting step towards 
constructing an efficient (in system size) and certified algo- 
rithm for local systems was recently made in Ref. ifTTl (com- 
pare Appendix G). 

Conclusions. A set of sufficient conditions for thermal- 
ization in quantum mechanics has been presented. The con- 
ditions are a natural translation of the standard assumptions 
from classical statistical physics. Along the way, a perturba- 
tion argument for realistic weak coupling has been proven that 
we expect to have significant applications beyond the scope of 
this Letter. By using our technical results we are able to de- 
sign a quantum algorithm preparing thermal states of arbitrary 
Hamiltonians with rigorous runtime and error bounds. 
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This is achieved using r ancilla qubits and running the algo- 
rithm an average number of 



Uuns < 2«J— e^+'^l!^-''" 



^ + B (13) 



times, where each run requires the application of n + 2r 
Hadamard gates, r controlled single qubit gates, n + q ( with 
q < r) single qubit measurements, and 2^ controlled unitary 
time evolutions under Hq = Hs + Hb for a time l/||i/o|loo- 

Notice that the time evolution under Hb can be imple- 
mented with m gates as the bath is a model of uncoupled 
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A Reduced density matrix of the dephased rectangular 
decoupled state - tlie general case 

In this section the distance between the subsystem of the 
dephased rectangular decoupled state and the Gibbs state 
'Pcibbs) computed. The discussion is similar to 
that in Ref. {6\. As all eigenvectors of Hq are given by 
tensor products of the eigenstates of Hs and Hb, \En) — 
\E^) (g) \E^), it is possible to trace the degrees of freedom of 
the bath of the decoupled rectangular state to get 



k \ 1 



(14) 



k=l 



with 



Pk 



EtU^liE~Ei) 



(15) 



Here, fl^{E — E^) denotes the number of states of the bath 
with an energy contained in the interval [E~E^, E—E^ +^]- 
For simplicity we assume that ||i?5||oo < E. Notice that the 
probabilities {pk} only depend on the shape of the spectrum 
of the bath. We aim at bounding the trace distance 

ds 

fc=l 



where 



Qk 



(17) 



are the eigenvalues of the Gibbs state and Zq := 

Let us define a function 5* as the logarithm of the number 
of states of the bath 



S{E) logng{E). 



(18) 



The Hilbert space of the bath is finite dimensional and thus the 
spectrum of the bath is discrete, from now on it is however as- 
sumed that the bath is large enough such that S : M"*" — M"*" 
can be well approximated by a twice differentiable function 
s : M+ — > K+. For natural systems and reasonable energy 
ranges the additional error from this approximation is expo- 
nentially small in the size of the bath. For the sake of concise- 
ness we ignore such exponentially small errors. We discuss 
this continuous approximation more closely in Appendix B. 

From now on we will hence focus on s{E) := logS£(_E), 
where S£ it a twice differentiable approximation of f2£. 
Taylor's theorem ensures that for every k there exists some 
^k^[E-E^,E], such that 



■<E-E^) 



<E)-^{E)El 



1 d'^s 
2dEP 



i^k)El 



= s{E) - f3Ei 



Ik: 



where 



is the inverse temperature and 



1 



Ik 



2dE^ 



i^k)El 



(19) 



(20) 



(21) 



A linear expansion of s{E — E^) is equivalent to an ex- 
ponential fit of the smoothed number of states of the bath 

^1{E-El), 



Thus, the probabilities {pk] can be written as 



Pk 



Eti^liE-Ef) 



(22) 



(23) 
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where Zp = ^ /3£^f+7i Therefore, the distance be- 

tween Lof^^^ and the Gibbs state depends on how well the 
density of states of the bath can be approximated by an expo- 
nential curve. The difference between the probabiUties reads 




(24) 



where the fraction Zq/Zp can be rewritten as 



Zq 
Zp 



k=l 

Introducing the notation 



~^P 



-Ik 



E 

fe=i 



-Ik 



7min := min min Iki^k), 
7max := max niax Jk{S.k), 



(25) 

(26) 
(27) 



we can write 



\Pk - qk\ 



Qk 



< qk (e^---^— -1) 



(28) 



The trace distance between the reduced dephased decoupled 
rectangular state Wn'"^'' ^nd the Gibbs state is thus bounded 
from above by 



1 



(29) 



where C is exponentially small in the bath size. We explicitly 
bound 7inax ^ 7min for a specific bath of uncoupled spins in 
Appendix B. 



B Reduced density matrix of the dephased rectangular 
decoupled state - a specific model 

In order to have a more explicit expression for the right 
hand side of Eq. (|29]l let us consider a particular model for the 
bath. We will also use this bath for the algorithm presented in 
the last part of the article. 

We start with the natural choice of m non-interacting un- 
coupled spin-1/2 particles with energies and 77. The spec- 
trum of this model is discrete and the eigenvalues are integer 
multiples of 77 between and ||i?B||oo = rim. The system 
is highly degenerate and the number of states with energy 77 k 
follows a binomial distribution (™) . 

This degeneracy makes it impossible to find a sufficiently 
good smoothed approximation for the number of 

states n-^{E) such that even for intervals whose width A 
scales like d^'^ with < k < 1 the error in the approxima- 
tion of the local density — 5l£(£'))/c? goes to zero for 



large (Ib- That is, the smoothed approximation would cause 
additional errors that would not go down exponentially fast 
with the bath size. 

Therefore, we need to assume that the degeneracy of the 
levels are lifted by a suitable perturbation. As we are only 
concerned about the spectrum and not the eigenstates of the 
bath, essentially any perturbation of adequate strength, such 
as basically any arbitrary weak interaction, will be sufficient. 

A convenient way to perturb the model for a theoretical 
study is to replace the fixed local field strength by a normal 
distributed random field strength such that on average the en- 
ergy of the local excited state is still 77 above the local ground 
state. If the width of the normal distribution is chosen in a 
suitable way, we get, with overwhelmingly high probability, a 
number of states n^{E) in the interval [E, E + A] that can 
be well approximated by 



fE+A 

^1{E) - / g''{E')dE' 
Je 



in the aforementioned sense, where 



g^{E') := -r 



77 



2 

7r77T, 



1/2 



-2m(^ 



(30) 



(31) 



Our aim is to bound 7inax — 7min in Eq. ( [29| l for this model. 
To do this, we first consider the case A <C H^^bHoo in which 
we can approximate as follows 



77 



2 

7r7TT, 



1/2 



(32) 



We can easily compute the second derivative of the logarithm 
of the above expression and find 



as 



(33) 



If A is not much smaller than ||i/B||oo, s{E) deviates from 
the parabolic form. But, larger A only make the curvature of 
s{E) smaller. In the general case we thus have 



> ^—-^ > 



and since j„ 
at 



dE^ 

<7n 



< 2\\H, 



s\ 



(34) 



J {rfm) we arrive 



< 



1 



-1 



(35) 



where C is exponentially small in the bath size. In natural 
situations the right hand side of the above equation will be 
very small as ||iJ5||^ <C 7^^777. 



C Rectangular states 



The rectangular states are a quite artificial class of states. 
So what happens if instead of Wp'''^'' we take states with a dif- 
ferent energy distribution? Will we still get something close 
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to a Gibbs state as long as the width of the energy distribu- 
tion is not too small and the density of states is exponentially 
increasing? The exponential increase in the density of states 
ensures that most of the contribution to ujj-^^ comes from the 
upper edge of the interval E + A in the case of a rectangu- 
lar distribution. Thus, it can be expected that one will get 
a state that is close to a Gibbs state for any energy distribu- 
tion that has a region with a bounded variation followed by 
a sufficiently sharp cutoff to higher energies. Where by suf- 
ficiently sharp we mean that it must drop from a value much 
larger than 1/d to a value much smaller than 1/d in an en- 
ergy interval that is small compared to (d'^s/dE^)^^/^. The 
bounded variation is required as for certain systems that vio- 
late the so called eigenstate thermalization hypothesis (ETH) 
ifTSl details of the energy distribution can have a huge impact 
on certain properties of the state. That bounded variation of 
the energy distribution is required to guarantee thermalization 
in such systems can be seen for example in the model studied 
inRef. 18J. 



D Complexity of the algorithm at low temperatures 

In this section we discuss the average number of times the 
algorithm must be run before the energy window correspond- 
ing to the desired temperature is hit. 

As we are only interested in the number of repetitions we 
can assume that the algorithm runs "perfectly" in that it gen- 
erates exactly a rectangular state with fixed width A at some 
position E of the spectrum. Furthermore, we are only inter- 
ested in cases where the reduced state of the rectangular state 
can be guaranteed to be close to a Gibbs state. Thus, we may 
assume that the spectrum of the bath is dense enough such 
that the number of states of the bath in the interval [E, i? + A] 
can be well approximated by a continuous and twice differen- 
tiable function S£ {E) that increases exponentially to higher 
energies. 

The probability of getting a rectangular state at energy E is 
given by 



P{E) 



ds 

^Eg{E-Ei), (36) 



fe=i 



where d — dsds is the dimension of the Hilbert space of the 
total system. Using that the density of states of the bath can 
be locally approximated by an exponential, we obtain 



PiE) 



1 ds 
k=l 



(37) 



Then, the probability of getting energy E can be bounded 



from above by 



P{E) 



< 



1 

ds 
1 



ds 

fc=i 



< 



(38) 



with E2 — Ef the gap of the system. The number of times 
that the program must be run on average in order to get the 
rectangular state at the position E is thus lower bounded by 



1 



> 



1 



PiE) 



j_ 

ds 



_^^-f3{Ei-E^) 



(39) 



This last equation shows that it is exponentially hard to go to 
low temperatures. This is to be expected, as no structure of 
the Hamiltonian is being used and generic local Hamiltonian 
systems can presumably not be efficiently cooled arbitrarily 
close to their ground state, because this would efficiently solve 
QMA-hard problems on a quantum computer 1 19 | . Needless to 
say, for specific Hamiltonians, with some additional structure, 
an algorithm can well be more efficient. The coefficient of 
the exponential is related to the features of the spectrum at 
low energies. In the case in which there is a gap, the ground 
state could encode the solution of a satisfiability problem that 
is expected to be hard to solve. 



E Temperature and its error for the bath described in 
Appendix B 

In this section we calculate how the inverse temperature (3 
and its error 6/3 depend on the energy value obtained while 
running the algorithm and the number of qubits q of register 
R that are measured. 

After running the algorithm a value is obtained and a 
state close to the rectangular state at energy ii^ = ||oo2^'s* 
is prepared. The system is then thermalized at an inverse tem- 
perature 



0- 



dE 



(40) 



Inserting the formula for (E) for the bath described in Ap- 
pendix B (l30| we find 



V V2 
4/1 

V V2 



E 
rjm 

2-« 



s* 1 



1 E 

\Hb\\od 



(41) 
(42) 



The algorithm will be run repeatedly until the value cor- 
responding to the desired f3 is obtained. For the bath consid- 
ered here this value is 



2« 



ll^sll 



\\He 



4 



(43) 
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It follows from Eq. (|43]l that 



V \ \\Hb\\oo 



(44) 



To reach a given precision 6/3 in the temperature it is thus 
sufficient to choose 



l0g2 



\H\\ 



log2 



Sl3r] 



\\Hs\\ 
\\Hb\\ 



(45) 



F Error and average number of runs of the quantum algorithm 

In this section the error of the quantum algorithm is derived 
and a bound on the average number of repetitions that are nec- 
essary to reach the desired energy is given. The total error is 
bounded by the sum of two errors 



with a more dense spectrum is harder This increase of the er- 
ror can however be compensated by increasing the accuracy 
of the phase estimation by linearly increasing the number of 
qubits r in the ancilla register. Similarly, lower temperatures 
can be reached by a quadratic increase of the number of an- 
cilla qubits of order 0{/3'^\\Hs\\'^)- As the phase estimation 
requires a total number of 2^ controlled unitary operations, 
this comes at an exponential cost in terms of runtime. This 
is expected to be a problem common to all algorithms using 
phase estimation on many particle systems. 

Before we go into the derivation of 'D{ujQc,L^n), let us 
quickly give the average number of runs that are required to 
get a certain temperature for the bath discussed in Appendix B 
in terms of A. From Eqs. ( |36l ) and ([39]l and the explicit for- 
mula for the smoothed number of states given in Eq. ( |73] l it is 
easy to see that 



runs < A^*^ 



+/3|lHslloo+' 



(52) 



2?Kc, PGibbJ <2?(WQC, 4°') +^?(w^^"^PGibbs)• (46) 

The deviation of the reduced rectangular state from the Gibbs 
state 'D{uJn^^\ Pqi^^^) has been calculated already in Ap- 
pendix A and B 



2?(-n"^°\pSibbs)<^(e^^-l 



C. 



(47) 



Here we bound the deviation I?(cjqc, Wp''') of the final state 
of the algorithm from the rectangular state. Our final result 
for the bath discussed in Appendix B will be 

2?(wQc,Wn)<e^^+''"''-""=°+ « (48) 
X 2"-'+^ (1 + ln(2'-«)/7r2) + C 

(again C is exponentially small in the bath size). It is conve- 
nient to write rj in terms of a dimensionless parameter A > 



1/2 



\H 



Siloo ■ 



because then the two errors can be written in the form: 



T>{ujQc,ujn) < 



J+/3||Hs|U + 



(49) 

(50) 
(51) 



X 2"-'-+^ (1 + ln(2'-«)/7r2) + C 

To reduce the trace distance 2?(wn*^°'' , PQibbg) it is favorable 
to choose a large A. Intuitively this captures the fact that the 
energy content in the bath must be much larger than that of the 
system in order to get a Gibbs state. However, at the same time 
this increases the error in the approximation of the rectangu- 
lar state, because performing the phase estimation on a system 



Deviation from the rectangular state for a general bath 

Here we bound the trace distance 'D{ujQc,uJn) between the 
state ujQc generated by the circuit with r ancilla qubits when, 
after measuring q of them, the outcome is obtained, and 
the rectangular state Wn in the interval [E, E + A] with E = 
s*2-9||i/||oo and A = 2«-'^||i?||oo is prepared. 

This trace distance can be written in terms of the one norm 
distance of the distributions of the diagonal elements 



1 

V{ujQC,UJn) = T^^kk -Pk\ , 



(53) 



k=l 



where the {qk} and {pk} are the eigenvalues of ujqc and Wn 
respectively, which can be written as 



Qk 



Pk 



Fr.q{(pk - <f) 
Gqifk - 'P) 

Zg 



(54) 
(55) 



where ipk Ek/\\H\\oo and (p :— E/\\H\\oo are the phases 
corresponding to the eigenvalues of the Hamiltonian and the 
energy of the rectangular state, Zp -^r.qi^j — 'f) 

and Zg :~ ^qifj ~ 'f) normalization constants, 

and the functions Fr^q and Gq are defined as 



fc=0 



Gqi^) 2' (eM-e((^-2-*)) 



with 



.sin^(7r2» 



(56) 



(57) 



(58) 
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proportional to the dimension of the Hilbert space d or, at 
least, to some power d'^ with < k < 1 such that L scales 
exponentially with the size of the bath. 

As h is continuos differentiable except at E and E + A, 
Taylor's theorem ensures that for all bins j (except for the 
two bins that contain E and E + A, which we can simply 
ignore because they only contribute an error of order l/L) 
there exists some value G [j\\H\\oo/L, {j + l)\\H\\oo/L[ 
for which 



2-1 



Figure 3. Plot of the Gq and Fr,q functions. As r increases, Fr.q 
tends to the rectangular function Gq . 



and B the step function (see main text and Fig.[3]l. 
Eq. ( |53] l is bounded by 

1 

■D{ujQC,ujn) < ^ X! -(p)~ Gq{ipk - 



where we have used that 



Fr,q{x) Gq{x) 



(59) 

1 1 

Zp Zg 



and 



Zg\ < 



k=l 



\Fr,q{Vk -V)- Gq{ipk - | 



Notice that Zq = Oa(£')2'', where Qa{E) is the number of 
states of the total system in the interval [E, £' + A]. 

Using that the spectrum of the whole system is sufficiently 
dense it is possible to approximate the upper bound in Eq. (|59]l 
by an integral. To make this rigorous we decompose the sum 
$9\ into bins of width \\H\\ac/L, 



L-l 



^=0 E,e[i,^[ 



(60) 



where 



h{E') 



1 



F^ 



E' 



\H\\ 



Gn 



E' 



\H\\ 



(61) 

is a function introduced to simplify the notation and L is the 
number of bins in which the spectrum has been divided. The 
idea is to take an L as large as possible such that the number of 
energy values in a bin VL\\H\\ac. / L{j\\H\\oo / L) can still be well 
approximated by its smoothed version S 1 1 1 1 ^ / ^ ( j 1 1 1 1 oo /i ) 
(compare the discussion of the twice differentiable approxi- 
mation in Appendix A and B). For baths for which such a 
smooth approximation is possible, one can always take an L 



h{Ek) = h 



mw 

L 



J\\H\ 



L 



(62) 



where j\\H\\^/L < Ek < {j + l)||iJ||oo/i- Then, the con- 
tribution of the j-th bin in Eq. (|60| can be bounded from above 
by 



L 



(63) 



- sup h'{i^\\H\\ 



L 



The last term in the parenthesis in Eq. ( [63) decreases with L 
and is thus exponentially small in the bath size. This can be 
verified with a lengthy, but straightforward calculation. 
Notice that the number of states can be written as 



Oii^ii^ {E') = p(E') 



\H\\ 



L 



+ 0{L-^). 



where 



AE' 



(64) 



(65) 



Putting Eqs. (|60|, ( |63| l and ( |64| l together, the error of the 
circuit can be bounded from above by 



iWHW 



i=t) 

0(L-' 



i\\H\ 



L 



\H\ 



(66) 



Notice that the upper bound of the previous equation con- 
verges to an integral for L — )■ oo with deviations that scale 
as 0(L"^). Thus, the trace distance between loqc and ujn is 
bounded from above by 



I?(wQC,^n) < \\hg\\i + 0{L-^) 



where 



\hg\U = 



h{E')Q{E')AE' . 



(67) 



(68) 



We will bound || /i f?]! i for the model described in Appendix B 
in the next section. 
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Deviation from the rectangular state for the bath described in 
Appendix B 

Next, the error of the circuit for a bath described in Ap- 
pendix B is bounded. The one norm p|| i from Eq. ( |67] i can 
be bounded by 

\\hQ\\i<^^A\Fr,q^G,h sup q{E'), (69) 

"Al-tij 0<£;'<||H||oo 



where the definition of h given in Eq. ( [6T| l has been used. We 
now find an upper bound on the supremum and a lower bound 
on the smoothed number of states. 

The density of states of the total system can be written in 
terms of the density of states of the bath as 



ds 



e{E')^Y.3''iE'-ES,). 



(70) 



k=l 



Inserting the density of states of our bath ( (ST) in the previous 
expression and writing it as a function of the inverse tempera- 
ture given in Eq. (|42]l we get 



k=l 



X df 



1 f 2 



77 \7rTO 



1/2 



(71) 



As we only consider positive temperatures /? > it is easy 
to see that the density of states attains its maximum value at 
/3 = 0, thus 



sup e{E') <^( — 
o<E'<\\H\\^ V K'^'m 



1/2 



(72) 



where d = ds ds is the dimension of the Hilbert space of the 
total system. 

The smoothed number of states Sa(£') can be lower 
bounded as follows 



fE+A 

Sa(S)= / g{E')dE' > g{E)A 

J E 



(73) 



•q \TTm 

where again we have assumed that is in a position of the 
spectrum with positive temperature. 
Finally we obtain 



+ 0{L-') 



\H\\ 



(74) 



In order for this error to become small the one norm \\Fr^q — 
Gq\\i must be made small enough such that it compensates the 
exponential prefactors. As we will see in the next section this 
can be achieved by using a polynomially large ancilla register 
R. 



One norm between F, .„ and G„ 



In this section we bound the one norm between Fr,q and Gq 
introduced in Eqs. (|56| and ( [ST] ). The one norm between Fn 
and G is defined as 



\H\\ 



F 



r,q 



E' 



— G„ 



E' 



\H\\qo J Vll-H^lloo 
By a simple change of variables it is easy to show that 

\\Fr,q - GqWl _ '""^ 



dE'. 



\H\\ 



\Fr,q{^)-Gq{^)\d^, 



(75) 



Remember that both ^ and Gq are normalized on the in- 
terval [0,1] and that Gq{ip) is a step function that is non- 
zero only in the interval [0, 2^''[, and that in this interval 
Fr.qif) < Gq{(p) — 2'^. Using this, the previous integral 
can be rewritten as 



\Fr_ 



Gq\\l 



\H\\ 



FN{ip)d(p 



(76) 



2""^' 1 



dip . 



Due to the symmetry and periodicity of f^, the contribution 
to the previous integral of the right tail of fr{ip — k2~^) is 
the same as the contribution of the left tail of fr{ip — (2^'' — 
fc2-'')) for/c2-'' < 2-9, thus 



fc2-'- + i 



2-9 
1 



fr{^-k 2-'-) dip 



(77) 



This implies that 



fr(p-2-i + k2-'')dp. 



\Fr,q — Gq\\l 



= 2"- 



E 

fe=0 



2-1 



fr {(p-k 2-'') dp . 



(78) 



The integral of the previous equation can be bounded by 



fe2-''+i 



where the bound 



fr{p-k2-')dp< 



cot(7r(2-9 - /c2-'')) 



7r2'- 



fr{p-k2--) < - 



sin-^ (p - A: 2-'') 



(79) 



(80) 



has been used. Therefore, the distance between the functions 
g and Gq can be bounded by 



\Fr,q Gq\\i ^ 2q-r+2 
lli/lloo 



1 ^ ^"^^ ^ cot (7r(2-9 - k 2-'^)) 



fe=0 



7r2'' 



(81) 
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where the sum has been split up in the k = 2^~'' case, which 
is exactly 1/2, and the rest. The sum of the previous equation 
can be bounded in two steps by 



2'-" 



^7r2-''cot(7rA:2-^) 



k=l 

< 7r2-''cot(7r2-'') + / cot(7ru)du (82) 

h 

<l + ln(2'-«). (83) 
Thus, the one norm between F/v and G is bounded by 

"^--^"^<2— + 4 + (84) 



(85) 



Moo \2 TT^ 

Inserting this into Eq. ( |74] i, we finally get 



I?(wqc, wn) < 



X 2«-''+2 (1 + ln(2'''-«) /tt^) + C , 

where C is exponentially small in the bath size. This com- 
pletes the discussion of an upper bound of the error made in 
the quantum algorithm. 



This is not made explicit in the paper, but appears to be cru- 
cial. Intuitively speaking, if a becomes larger, the dephasing 
is more exact, but then (6)-(8) can no longer be used. 

This step is then used in the procedure following Eq. (13). 
A C is introduced and dephasing between eigenstates with rel- 
ative gap larger than ( is considered. Then a new Hamiltonian 
H is constructed, which has the following feature: It has the 
same eigenbasis as iJ + eh, but with eigenvalues grouped in 
bins, such that the smallest gap between bins is (. For the 
following procedure to work, one has to take (again in the 
asymptotic notation) 



a = 6(1/0, 



(87) 



so both (7 = 0(1/C) and a = ^{1/0 (although only a = 
0{l/() is made explicit, both bounds are actually needed). 

This, however, seems to already fix ^ = ri(e). So there is 
no longer the freedom to have 



(88) 



which would mean that ( = ^(e^)- This appears to contradict 
the above statement. Again, intuitively, a is forced to be small 
for the Dyson approach to work, but at the same time this 
gives a constraint to C, which requires a to be large. 



G Discussion of the argument presented in Ref. flTl 

Ref. ifm presents a novel approach towards thermalizing 
quantum systems using an iterative approach, in which pre- 
thermalized parts are put together in a suitable fashion in or- 
der to arrive at a Gibbs state of a quantum system with a local 
Hamiltonian. This argument provides a new intuition on how 
one can think of thermalizing local quantum systems, differ- 
ent from a quantum Monte Carlo approach (and the one pre- 
sented here). In this appendix, however, we point out a subtle 
mistake of the published version of the paper; however, this 
challenge can be overcome and a new proof can be formu- 
lated that retains key elements of the previous structure of the 
argument |23|. The intuition behind the argument and in fact 
the algorithm itself hence turn out to be correct. 

Each merging step consists of two steps. The first is a prob- 
abilistic step that updates the probability weights of the Gibbs 
state by means of postselection. The second one aims at rotat- 
ing the eigenbasis of the old Hamiltonian to the one of the new 
Hamiltonian by means of an instance of dephasing. Each step 
has as input a chosen e > 0. For the entire algorithm to work, 
this procedure has to be correct up to errors of O(e^). In what 
follows we refer to the equation numbering of the preprint v2. 

In Eq. (5), perfect dephasing is being achieved when a — > 
oo. This is approximated by imperfect dephasing with a finite 
(J. From the Dyson series of second order (6)-(8), it follows 
that (in the asymptotic notation) 

CT = 0(l/e). (86) 



H Tlieorem 1 for an exponential density of states 

As we have seen in the main text in theorem 1 the distin- 
guishability of the microcanonical states lJ^^^ and Wn corre- 
sponding to an interval [i?, i? + A] of the Hamiltonians 
and R = Hq + V is bounded by 

2'(wn,Wn ) < T>{ujn,ujn ) < ^ -7^ , (89) 

where il,„ax and Aft are the maximum^ and the difference, 
of the dimensions of the support of lu^ and cun, and fl^ is 
the total number of eigenstates of H and Hq in the intervals 
[E, E + e]and[E + A - £,E + A]. 

In the main text, in order to give a more comprehensible 
interpretation to Eq. ([89|, an approximately constant density 
of states was assumed and it was shown that 

2'(-n^-n^^°^)^^(&^)'''. (90) 

However, we have seen that thermal states emerge in situa- 
tions where the density of states is locally well approximable 
by an exponential. Thus, new conditions that ensure the in- 
distinguishability of the microcanonical states a;^ and Wn'^"'' 
must be derived for the exponential density of states 

g{E) cx e^^ . (91) 

In order to do this, let us notice that both terms in the upper 
bound of Eq. ([89| are positive and must be simultaneously and 
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independently small. For the first term, this is only possible 
if 1 1^1 loo ^ £■ To find the condition for the second term, let 
us assume that the interaction does not shift excessively the 
energy levels such that AO/ilmax can be neglected. Then, 



1 > 



> 



2(1 



3-/3A 



)' 



(92) 



and the condition (3e 1 is required. These two necessary 
conditions can be summarized as 



(93) 



It is easy to see that they are also essentially sufficient condi- 
tions by using 



< 



rE+A 

Ie+a- 



, giE')dE' 



1 - c 



j-E+A 

Je 



g{E')dE' 



-PA 



< 



(3e 



-/3A ' 



to bound Eq. ((89| as follows 



S(0)^ 



< 



Af2 



1 - e"/3A 



(94) 
(95) 



Finally, let us simply choose e = \J\\V\\oo/ P- The trace 
distance between the microcanonical states Wp"^ and Wn then 
reads 



< 



1 



-.-IIA 



(96) 



Equation ( [96[ i ensures the indistinguishability of the inter- 
acting and non-interacting microcanonical states as long as 
A > fc^T is not too small and the condition 



(97) 



where ks is the Boltzmann constant and T the absolute tem- 
perature, is fulfilled. That is, Eq. ( |96] l gives us a physical 
intuition about when an interaction is weak in the sense of 
theorem 1 : An interaction is weak if it is small compared to 
the thermal energy ksT, which is a measure of the intensive 
energy content of the system. In conclusion we see that how 
strong an interaction feels for a system depends on how much 
energy it contains. 



